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ABSTRACT 

We present an analytic solution to diffusion equation for high energy cosmic 
rays in the expanding universe. The particles are assumed to be ultra-relativistic 
and they can have energy losses arbitrarily dependent on energy and time. The 
obtained solution generalizes the Syrovatsky solution, valid for the case when 
energy losses and diffusion coefficient are time-independent. 

Subject headings: extragalactic cosmic rays, diffusive propagation of cosmic rays. 


1. Introduction 

The cosmic rays is important component of extragalactic medium, which is responsible 
for heating of extragalactic gas and for production of various types of radiation in the space. 
The charged cosmic ray particles (electrons, protons and nuclei) propagate diffusively in ex¬ 
tragalactic space due to scattering in the magnetic fields. The regime of diffusive propagation 
is reached if scattering length, £ sc , is much less than the size of the considered region. 

Typically, diffusion takes place in the magnetized plasma, where the particle deflections 
occur due to scattering off the turbulent pulsation and hydromagnetic waves. This process 
proceeds in the resonant regime (see (Lifshitz & Pitaevskii 2001), Sections 55 and 61), when 
a particle giro-radius is equal to a wave length. The calculation of the diffusion coefficient 
for the resonant scattering of particles off the hydromagnetic waves is presented in the book 
by Berezinsky et al. (1990a). 

The value of diffusion coefficient is quite different for three distinct regions in extragalac¬ 
tic space: voids, filaments and clusters of galaxies. The diffusion in these regions occurs only 
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at energies when the diffusion length, I h iff , is much smaller than the size of the region L. In 
the opposite extreme case particles propagate (quasi)rectilinearly. 

For propagation of high energy particles from a single source at point r g the diffusion 
equation reads 

§i n p( E i 'F *) “ div [ D ( E , r, t)Vn p ] - [b(E , r, t)n p ] = Q(E, r, t)S 3 (f - f g ), (1) 

where n p (E, r, t ) is the space density of particles p with energy E at time t and at the point 
r, D(E,r,t) is the diffusion coefficient, b(E,r,t) = —dE/dt describes the continuous energy 
losses, and Q(E,f,t ) is the source generation function. 

In the case when D, b and Q depend only on energy, the method of exact analytic 
solution to the diffusion equation has been suggested by Syrovatsky (Syrovatskii 1959). For 
a single-source diffusion equation (1) the solution for spherically-symmetric case can be 
presented (Berezinsky et al. 1990b) as 


n p(E,r) 


1 

W) 



dE g Q(E g 


exp 


4A (E,E a ) 


[4ttA (E,E g )} 


3/2 


( 2 ) 


where 

KE,E.) = j '*™| (3) 

is the Syrovatsky variable which has the meaning of the squared distance traversed by a 
particle in the observer direction, while its energy diminishes from E g to E. 

The other Syrovatsky variable, 



( 4 ) 


has a meaning of time, during which the particle energy diminishes from E g to E. 

In this paper we shall present an analytic solution to the diffusion equation in the 
expanding universe for the case when D, b and Q are arbitrary functions of energy and 
time. We use the computation method which differs from that of (Syrovatskii 1959). While 
Syrovatsky used method of the Green functions, we solve the equation for a single source. Its 
position is described by the delta function <5 3 (r — r g ), and thus the equation for the Fourier 
component does not contain 5-functions. This is the first order linear equation with partial 
derivatives, which can be solved by a standard method, introducing auxiliary characteristic 
equation. 
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As far as physical applications are concerned, the solution for a single source is most 
general, because all other cases, e.g. with homogeneously distributed sources or with single 
non-stationary source (see (Berezinsky et al. 1990b)) can be straightforwardly obtained from 
it. 


2. Diffusion equation in expanding universe 


We shall use the Friedmann-Robertson-Walker metric for the flat space and radial di¬ 
rection, following Weinberg (1972) 

ds 2 = c 2 dt 2 — a 2 (t)dx = —g^dx^, (5) 


where diag g, JV = (—1, a 2 , a 2 , a 2 ) and diag g lw = (—1,1/a 2 ,1/a 2 ,1/a 2 ), x is the spatial 
coordinate, corresponding to comoving distance, and a(t) is the scaling factor of expanding 
universe, normalized as a(t o) = 1 at present age of the universe to- The redshift z is given 
by 1 + z = 1 /a(t) and dt/dz by 

dt 1 

dz H 0 ( 1 + z)yfVL m (l + z ) 3 + A 

where H 0 is the Hubble parameter at z — 0 and Q m and A are cosmological mass density 
and vacuum energy in units of the critical density. 

The physical and proper distances are locally determined as dr = a(t)dx. For the 
proper distance it is assumed that a(t) is not changed in the process of distance measurement 
between x = 0 and x, and thus the proper distance between these two coordinates is r prop = 
a(t)x , and velocity of the universe expansion is u = cix = H(t)f prop . 


For the physical distance it is assumed that measurement is performed with help of the 
light signal (ds 2 = c 2 dt 2 — a 2 (t)dx 2 = 0) and the distance between the object with redshift 
z and observer with z = 0 is 


r p h 



a(t)dx = c 



dz 


dt 

dz 



dz 

H 0 (l + z)\J fl m (l + ^) 3 + A 


( 7 ) 


Following Peebles (1980) we shall use in this paper the proper distance, denoting it by 
r and r without subscript. Being the formal quantity, the proper distance coincides locally 
with the physical distance, and in most applications we shall in fact use only local properties 
of the proper distance. The positions of the particles will be described by both proper 
distances r and by comoving distances x, the latter can be considered as the set of fixed 
coordinates in the expanding universe. 
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Expansion of the universe affects diffusion, and we shall take it into account using the 
proper formalism. 

The local observer sees the particle flux density produced by diffusion 

= (8) 

In case of isotropic diffusion (the diffusion coefficient is rotation invariant), jk is the covariant 
space vector, which together with proper space density of particles n forms the covariant 
4-vector j^(x,t) = ( n,jk ) (the Latin indices run through 1-3, and the Greek ones through 
0 - 3). 

The conservation of current can be written as (Weinberg 1972) 

£: (-Tsf) = o. (9) 

where g = \det g lw and ^fg = a 3 (f). Differentiating Eq. (9) with using definition of the 
Hubble parameter H(t ) = a(t)/a(t), then transforming the contravariant component j k into 
covariant one as j k = g km j m , with diag g km = l/a 2 (f), and finally using Eq. (8) for 

A/ = (io) 


one obtains the diffusion equation 


t) + 3 H(t)n(x, t) - ^V 2 x n(x, t ) = 0, (11) 

where 3 H(t)n term describes expansion of the universe. 

The derived equation can be obtained from conservation equation (9.15) of Peebles 
(1980) excluding there the advection term and adding the diffusion term. 

The derivation of the diffusion equation above is rather formal, and now we shall present 
it in the transparent and physically clear way. 

We shall operate here and everywhere below with particle density n in the expanding 
volume using two sets of the coordinates (r,t) and (x,t). This density is given by 

n(r,t) — n(x,t), (12) 


with f(t) = a{t)x. Differentiating n( f, t), 


Mr, t) = 


dn \ dn dr / dn 
dt ) „ dr dt V dt 


+ H{t)r 
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and using dn(r,t)/dt = dn(x,t)/dt from (12) one obtains the Peebles relation given by 
Eq. (9.13) from (Peebles 1980): 


dn 

~dt 


dn 

~dt 


+ H(t)f'V r n, 


(14) 


where subscripts x and r indicate basis (x,t) and respectively. 

We shall obtain the diffusion equation from conservation of number of particles. Con¬ 
sider a sphere of radius x, which expands in the basis (r,t) as r(t) = a(t)x. The number 
of particles inside this sphere is changing only due to diffusive flux, which is defined as 
j = —D(f,t)'V r n(r,t). The corresponding equation reads: 

— f dVn(f,t) — — f jds= j dVdiv[D{f,t)V r n(r,t)\, (15) 

dt Jv(t) Js(t) Jv(t) 


where S(t) is the expanding sphere and V(t) is the volume inside it. In Eq. (15) the Gauss 
theorem was used. 


Performing differentiation in the lhs of Eq. (15) with respect to time, and taking into 
account expanding of the elemental volume dV with time, 

jW = j t [a 3 (f)(TP) comov ] = 3 H(t)8V, 

and using for dn/dt Eq. (13), one obtains 


dn 

~dt 


+ H(t)f'V r n + 3 H(t)n — div D(r, t) V r n(r, t) 


= 0 . 


(16) 


This is the diffusion equation in the physical basis (r, t). One may add to the rhs the source 
term Qo5 3 (r — f g ). 

The sum of the second and third terms in Eq. (16) merges into V r (nu ), where in our 
case u = H(t)r is expansion velocity of the universe, and in absence of the diffusion term 
we obtain Eq. (9.11) from (Peebles 1980). 

In the basis (x,t) the first two terms in Eq. (16) give dn(x,t)/dt (see Eq. 14), and we 
arrive at the diffusion equation in the form 


dn 

~dt 


+ 3 H(t)n(x, t) 


Rill\7 2 n(x t) = d 3 ( 

a*(ty x{,) a 3 (t) 1 


x 


9)i 


(17) 


in accordance with Eq. (11). Deriving Eq. (17) we assumed that D is a function of time 
only. 
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Eqs. (11) and (17) are the simplified diffusion equations without energy loss and full 
source terms, which should be included additionally. Introducing these terms as nb ) 

and Q(E, t)5 3 (f— f* g ) we arrive at the diffusion equation for expanding universe in the basis 
(x, t ) convenient for its solution: 


dn 

~di 


b(E, t) 


dn 

dE 


+ 3 H(t)n — n 


db(E,t ) 
dE 


D(E,t) 
a 2 (t ) 


V 2 r n 


Q(E,t) 

a 3 (t ) 


(S 3 (X - Xg), 


(18) 


where the total energy losses dE/dt = — b(E,t ) may be presented as a sum of collisional 
energy losses b in t(E,t ) and adiabatic energy losses H(t)E, and Q(E,t) is the number of 
particles with energy E produced at time t per unit time. In Eq. (18) b(E,t) and D(E,t) 
are the arbitrary functions of time and energy, and n = n(t,x,p) must depend in fact on 

Note, that n(t, x,p) in Eq. (18) is not the distribution function f(t, x, p ), which is defined 
in statistical physics (see e.g. (Lifshitz & Pitaevskii 2001)) as the density in phase space. 
The relation between them is given by n(t,x,p) = 47r p 2 f(t,x,p). We remind that here is 
considered the ultra-relativistic case p ~ E/c. 


3. Analytic solution to the diffusion equation 

We shall introduce the Fourier transformation 

n(t, x, E) = W J dufUE, t)e^~ 3 ’K (19) 

Using the Fourier expansion of the 5 function 

5 3 {x-x g ) = J due (20) 

we obtain from Eq. (18) equation for the Fourier components f w (E,t): 

f (F f) — ^ foi \ 

The characteristic equation for Eq. (21) is 


j t UE,t)-b[E,t)-Fj-JE,t) + 




dE 


a?(t) 


dE_ 

dt 


-b(E,t), 


( 22 ) 


with the solution S' = E'(E,t,t') being identical with the generation-energy trajectory 
E g (E,t,t') used in (Berezinsky & Grigorieva 1988) and (Berezinsky et al. 2002), where E g 
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is the energy with which a particle must be generated at time t' in order to have at t the 
observed energy E. Here and everywhere below we shall use the energy E at the time t of 
observation to mark a characteristic trajectory S' = E'(E,t,t'). Sometimes for brevity we 
shall omit t. 

The solution to Eq. (21) with energies taken on characteristics is given by 


fu(E,t) = j dt'^^y-ex p<j - J df 

in l tj 




d£ h 


a?(t 


(23) 


where t g is the generation time. 

We introduce now A (E,t,t'), the analogue of the Syrovatsky variable given by Eq. (3) 


a (*, o = !‘df D(£ "' th 


where S" = E"(E,t,t") is the characteristic trajectory. 

The exponent in Eq. (23) can be simplified using the notation 


0^(77, f, t') = / dt" 

Jt’ 


3 H(t" 


db(£", t" 
d£" 


+ cn 2 A (E,t') 


= In 


i±dj 

1 + z J Jt' 


d£' 


Thus, 


e ~a w (E,t,t') _ / \ + ~ \ e -<l5 2 A (E,t r ) 


[' , ,,db(£\t n 
exp / dt” v ’ 

Jt' 


1 + z’J " " ()£" 

Coming back from the Fourier component f w to density n, using identity 


iu(x — x g ) — ci; 2 A = —A 
and performing integration 

[ du exp 




LU — l- 


2A 


(X - Xn 


X . - \X-Xg) 

-A w - i- 


2A 


4A 


= (*/A ) 3/ t 


we obtain the solution as 


n(t, x, E) = 


ir 


3/2 


(2ir) 


ex P 


0 X-Xg ) 2 

' 4A(E,t') 


(1 + ^) 3 [A (E,f)] 3/2 


exp 




tv 


B£‘ 


(24) 


(25) 


(26) 


(27) 


(28) 



Now we shall find the solution for t = to (2 = 0 ) changing the variables t' —> z, t" —> z! and 
presenting the energy loss as the sum of adiabatic and collisional (interaction) energy losses: 

b(E, t ) = H(t)E + b int (E, t ). (29) 


It gives 


n(t 0 ,x,E) = 


[ Zg 

dt 

/ dz 

Jo 

dz 


exp[ 


(1 + z)Q(S g , z) exp 


(x - x g ) 2 /AX(E,z)\ 


-1 

is* 

dt' 

db int (£', z'Y 

1 U//C- 

L^o 

dz' 

d£' 


x 


[4ttA (E,z)] 3 / 2 


(30) 


where £ g = E g (E, z) is the generation energy in the source taken on the characteristic. 

In this expression one can easily distinguish the term dE g /dE , given for the case of 
protons interacting with CMB in (Berezinsky & Grigorieva 1988) and (Berezinsky et al. 
2002). In our case it equals to 


dE g 

dE 


(1 + z) exp 



db int (£',z')\ 
d£' ) ' 


(31) 


Now the solution for the case of arbitrary energy losses b int (E,z ) can be written in the 
compact form 


n(t 0 ,x, E) 



Q[E g (E,z),z 


exp [—(x — x g ) 2 /AX(E, z)] dE g 
[4ttA (E,z)} 3 / 2 lE’ 


(32) 


where n(t 0 ,x,E ) actually depends on x — x g as was foreseen above. 


In the case of ultra-high energy protons interacting with CMB, dbi nt jdE is given ac¬ 
cording to Berezinsky & Grigorieva (1988) and Berezinsky et al. (2002) by 


dbj nt (E, t) 
dE 


= (1 + ^) 2 


db 0 (E') 
dE' 


(33) 


E'=(l+z)Eg{E,z) 


and 


dEg_ 

dE 


= (1 + z) exp 


dz 


(1 + *) 5 


Ho J 0 1 + z ) 3 + A 


db 0 (E') 
dE' 


E'={l+z)E g {E,z) 


(34) 


where b 0 (E) is the interaction energy loss at z = 0. 

Eq. (32) presents our final result in the form of the particle density n(to,x,E) in terms 
of comoving distances x for arbitrary energy loss bi n t(E,t ) and diffusion coefficient D(E,t). 
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The main difference with the Syrovatsky solution is due to A (E,t) given by Eq. (24) and 
dEg/dE given by Eq. (31). Lemoine (2005) has heuristically written a solution to the dif¬ 
fusion equation in expanding universe, including only adiabatic energy losses. He used the 
Syrovatsky solution for the diffusion equation in terms of conformal time rj. However, the 
Syrovatsky solution is not valid there, because both diffusion coefficient and the Hubble pa¬ 
rameter depend on conformal time rj. Besides, our equation (18) differs from the Lemoine’s 
equation by extra a(t) in the diffusion term. 


Eq. (32) becomes the Syrovatsky solution (2), when D(E,t) and b(E,t) do not depend 
on time and a(t) = 1. In this case dE g /dE = b(E g )/b(E) (Berezinsky & Grigorieva 1988), 
and using b(E g )dt = dE g one obtains Eq. (2) with 




being the Syrovatsky variable (3). Proper transition to the Syrovatsky solution is one of the 
tests of our solution (32). 


4. Tests of the method 

We want to test our method of equation solution for some cases, when solutions are 
known. One of them is given by the case of time-independent energy loss and diffusion 
coefficient, which must result for a(t) = 1 in the Syrovatsky solution. In Section 3 it has 
been already demonstrated that our solution passes this test. 

The second test consists in the convergence to the universal spectrum, and the third, 
most important one, examines the cosmological part of our solution: The solution of equation 
for the rectilinear propagation is obtained by the same method as for diffusion and result 
coincides with well known formula. Below we shall describe these tests. 


4.1. Convergence to the universal spectrum 

According to the propagation theorem (Aloisio & Berezinsky 2004) in the case when 
distances between sources (d) become smaller than propagation and interaction lengths (e.g. 
when d —> 0), the diffuse spectrum has an universal form, independent of the mode of 
propagation. 

For the power-law generation spectrum with exponent r ) g > 2 , 

Q(E,z) = (7 S -2)L„(1+2)”E-V 


(35) 
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where Lq is the particle luminosity of a source at z = 0 and a takes into account a hypo¬ 
thetical cosmological evolution (in Eq. (35) all energies are measured in GeV, luminosity in 
GeV s -1 and E min = 1 GeV). In this case the universal spectrum for the diffuse flux J(E) is 
given (Aloisio & Berezinsky 2004) as 


J(E) 





dt 

dz 


(1 + z) m Ef°{E,z) 


dEg 
dE ’ 


(36) 


where £ 0 = L 0 n s (Q ) is the emissivity at z = 0, n s (0) is density of the sources at z — 0, 
n s (z) = n s (0)(l + z) 13 describes the hypothetical source evolution and a + (3 — m. 

Let us now calculate the diffuse flux from the density of the particles n(to,x,E) given 
by our general solution (32). According to the propagation theorem we must obtain the 
universal spectrum (36). 

Assuming the homogeneous distribution of the sources n s (z) in the space (which provides 
the propagation theorem) we can find the diffuse flux multiplying Q(E g ,z ) in Eq. (32) by 
the density of the sources n s (z) = n s ( 1 + z) 13 and integrating over positions of the sources in 
the coordinate space. Using (35) we have 


J P (E) = ~ 2 )A) / ^ x ldxs / dz 


dt 


dz 


(1 + z) m Ep»(E,z) 


exp[— x 2 s /A\{E, z)\ dE g 


[4vrA(£, z)f/ 2 dE ’ 
(37) 


where x s = x — x„. Changing the order of integration and using 


dx 


inx 2 


o (4vrA) 3 /2 


exp - 


x 

4A 


= 1 , 


we arrive indeed at the universal spectrum (36), as it must be. The assumption of the 
power-law generation spectrum in Eq. (35) does not reduce the generality of the proof. 


4.2. Rectilinear propagation 

The third test which we shall study here is less trivial than other two and it examines the 
cosmological part of our solution. We consider the rectilinear propagation of ultra-relativistic 
particles. This case is well known, and is given by light propagation. We shall find a solution 
for this case, solving the propagation equation by the method employed in Section 3. 

The rectilinear propagation of particles with velocity v results in adding this velocity 
to the expansion velocity H(t)r in the rhs of Eq. (13) and in excluding the diffusion term 
from the equation. Introducing the unit vector e in the direction of propagation, and using 
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dr = a(t)dx we obtain from Eqs. (13) and (18): 


dn ce dn dn db Q(E,t).n,_ _ 

— —I—— b(E, t) —— + 3 H(t)n — n—— = — S 3 (x — x„), 
dt a(t) dx y ’ ’ dE w dE a 3 (t) K 9h 

where n = n(t, x, E). 

Performing the Fourier transformations (19) and (20), and using 


(38) 


ce d . _ . 1 

; n(t, x, E) = 


du ^ C f u (E,t)e 
a(t) 


iu(x—x g ) 


a(t ) dx v ’ ’ ' (2vr) 3 

we obtain from (38) the following equation for the Fourier transform f w (E,t): 




3 H(t) - 


dE 


a{t) 


U e A = ^M L - ( 39 ) 

a 6 [t) 


With the help of the characteristic equation (22), which gives S' = E'(E,t,t'), the solution 
of Eq. (39) is found as 

/„(£,«)= «,<■)_ 

Jtn a ) 


with 


a w (E,t,t') = / dt" 


3 H(t" 


db(S",t " 
dS" 


+ i 


ecu c 


a(t" 


(40) 

(41) 


To calculate cc aJ (F,t,P) we use 


rt i + z ' 

dt"H(t") = ln^-, 

4 1 + Z 


^ dt" 


z ecu c 


a(t" 


= iu[x(t) — x(t')], 


and b(E,t ) = H(t)E + bi nt (E,t). Then it follows 

2 




1 + z 


1 + z' 

Coming back to n(t, x, E) we have 


e -H3{x-x') exp 


1 

M 

dt" 

1 

dz" 


db int (S", z" 
dS" 


(42) 


n(t, x, E) = 


du 

( 2*0 


_ e «3(3-3 a ) ( 1+z X 


a 3 (t') \1 + P 


x 


exp 


I dt" 
dz" ' ^ 


dz" 


db int (S",z " 
dS" 


(43) 
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Using the Fourier expansion of S function and its properties, 

f J^_ e ^- 3g ) = <$3^) _ = _l^S[x(t') ~ X g \ 


J (2 7T 

we get the solution at t — to (z=0) and x = 0 
n(to, x — 0, E) — 


4 ttx 2 


a (tg) 

4:7lCXg 




dt' Q(£ ': t!) a ^ 9 \ d(t' — t g ) exp 


1 

dt" 

dbint(S", z") 

1 U//C- 

1 Jo 

dz" 

dE" 


a(t') Aircx 2 
where a(t')( 1 + z') — 1 is used. 

Performing integration over f and using Eq. (31) we finally have 


(44) 


n(t 0 , E) = 


Q(Eg,tg ) dEg 

47rca;g(l + z g ) dE ’ 


(45) 


where x g is the comoving distance to a source. 


One can compare n(to, E) from Eq. (45) with energy flux F (in erg cm -2 s -1 ) of photons 
emitted by a source with luminosity L at comoving distance x g , as given by Eq. (2.42) in 
the book by Kolb & Turner (1990): 


F = 


(46) 


4ra|(l +%) 2 

where one factor (1 +z g ) arises from the time dilation and the other one from energy redshift. 

In our case L —> Q(E g )dE g and F —> cn(E)dE, with one factor (1 + z g ) disappearing 
because we consider the number of particles instead of the luminosity and energy flux. 


Therefore, our method passes this test, too. 


5. Conclusions 


We shall conclude giving the formulae in the form convenient for practical use. 


Our basic diffusion equation (18) for ultra-relativistic particles propagating from a single 
source, reads 


dn 




+ 3 H(t)n — n 


db(E,t) 

dE 


D(E,t) 
a 2 (t ) 



Q(E,t) 
a 3 (t ) 


5 3 (x-x a ), 


(47) 


where n(t, x, E) is the particle number density per unit energy in expanding volume of the 
universe, x is coordinate corresponding to the comoving distance, dE/dt = —b(E, t ) describes 
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the total energy losses, which include adiabatic H(t)E and interaction bi nt (E , t) energy losses, 
and Q(E , t) is the generation function, given by the number of particles generated by a single 
source, located at coordinate x g , per unit energy and unit time. 


Solution of Eq. (47) can be presented in the spherically-symmetric case as 

exp[-Xg/4A(£, z)] dEg 


n 


(Xg, E) = 


r z a 

dt 

/ dz 

Jo 

dz 


Q[E g (E, z), z 


where 


dE 


A {E,z) = 


= (1 + z) exp 


[ dz' 

dt' 

Jo 

dz' 


[An\(E,z)} 3 / 2 

D(£',z') 

a 2 (z') 


dE 


-1 

dt' 

db int (S\z'Y 

1 Uj/J 

Uo 

dz' 

dS' 


(48) 

(49) 

(50) 


where S' = E'(E,z') is a characteristic trajectory, which gives energy E' of a particle at 
epoch z', if this energy is E at z = 0; E g (E, z) has the same meaning. The upper limit z g in 
the integral of Eq. (48) is provided by maximum energy of acceleration as E g (E, z g ) = E max , 
or by z max , what is smaller. 


The solution (48) is intentionally presented in the form similar to the Syrovatsky solution 
with A (E, z) being an analogue of the Syrovatsky variable. However, this solution cannot be 
obtained neither from the Syrovatsky solution (2) nor by the Syrovatsky method of solution, 
which is based on time-independent quantities A and r as new variables. 


Eq. (48) is convenient for various applications. 

In case the sources are distributed homogeneously with density n s (x g ) the total density 
of particles n tot (E) can be found by integration of Eq. (48) over Anx 2 g dx g , with the diffuse 
flux given by J(E) = (c/An)nt ot (E). In case of discrete distribution of the sources, the 
density n to t(E) must be found by summation over sources, like it is done in (Aloisio & 
Berezinsky 2004). For non-stationary source, being switched on at redshift z\ and switched 
off at redshift Z 2 , the limits of integration in Eq. (48) must be taken as z\ and Z 2 with 
additional acceleration cut on z max given by equation E g (E, z g ) < E max (see also Berezinsky 
et al. (1990b)). 
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